Interdependent networks: Reducing the coupling strength leads to a change from a 

first to second order percolation transition 
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We study a system composed from two interdependent networks A and B, where a fraction of the 
nodes in network A depends on the nodes of network B and a fraction of the nodes in network B 
depends on the nodes of network A. Due to the coupling between the networks when nodes in one 
network fail they cause dependent nodes in the other network to also fail. This invokes an iterative 
cascade of failures in both networks. When a critical fraction of nodes fail the iterative process 
results in a percolation phase transition that completely fragments both networks. We show both 
analytically and numerically that reducing the coupling between the networks leads to a change 
from a first order percolation phase transition to a second order percolation transition at a critical 
point. The scaling of the percolation order parameter near the critical point is characterized by the 
critical exponent (3 = 1. 
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Most of the research on networks has concentrated on 
the limited case of a single network [H-Q while real world 
systems are composed from many interdependent net- 
works that interact with one another [6|-l8j- As a real 
example , consider a power-network and an Internet com- 
munication network that are coupled together. The In- 
ternet nodes depend on the power stations for electricity 
while the power stations depend on the Internet for con- 
trol 0. 

We show that introducing interactions between net- 
works is analogous to introducing interactions among 
molecules in the ideal gas model. Interactions among 
molecules lead to the replacement of the ideal gas law 
by the Van der Waals equation that predicts a liquid-gas 
first order phase transitions line ending at a critical point 
characterized by a second order transition (FigQJa)). 
Similarly, interactions between networks give rise to a 
first order percolation phase transition line that changes 
to a second order transition, as the coupling strength be- 
tween the networks is reduced (Figfljb)). At the critical 
point the first order line merges with the second order 
line, near which the order parameter (the size of giant 
component) scales linearly with the distance to the crit- 
ical point, leading to the critical exponent (3 = 1. 

In interdependent networks, nodes from one network 
depend on nodes from another network. Consequently, 
when nodes from one network fail they cause nodes from 
another network to also fail. If the connections within 
each network are different, this may trigger a recursive 
process of a cascade of failures that can completely frag- 
ments both networks. Recently, Buldyrev et al [10j stud- 
ied the coupling between two N node networks A and B 
assuming the following restrictions: (i) Each and every 
node in network A depends on one node from network B 
and vice versa, (ii) If node A4 depends on node Bi then 
node Bi depends on node Aj. They show that for such a 
model when a critical fraction of the nodes in one network 
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FIG. 1: (a) The van der Waals phase diagram. Along 
the liquid-gas equilibrium line the order parameter (density) 
abruptly changes from a low value in the gas phase to a high 
value in the liquid phase. At the critical point(P c ,T c ) the 
order parameter changes continuously as function of temper- 
ature if the pressure is kept constant at the critical value, but 
its derivative (compressibility) diverges. This is a characteris- 
tic of the second order phase transition, (b) The percolation 
phase transition for two interdependent networks as obtained 
from the numerical solution of system ([7]) for qs — 1 and 
a — b — 3. Here 1 — p, the fraction of removed nodes from 
network A, plays the role of temperature. (As 1—p increases, 
the disorder increases.) The fraction 1 — qA of independent 
nodes in network A plays the role of pressure. (As 1 — qA in- 
creases the stability of network A increases.) Below the criti- 
cal point, the system undergoes a first order phase transition 
at which, (3oo, the fraction of nodes in the giant component 
of network B abruptly changes from a finite value to zero. As 
we approach the critical point, /3oo — ¥ 0. Above the critical 
point, the system undergoes a second order transition where 
the giant component continuously approaches zero. 



fail, the system undergoes a first order phase transition 
due to the recursive process of cascading failures. 

However, when examining the features of real inter- 
dependent networks such as the power network and the 
communication network presented above, we observe that 
in practice not all nodes of network A depend on network 
B and vice versa. We therefore introduce a general model 
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that is applicable to many real networks. The model con- 
sists of two networks A and B with the number of nodes 
N A and N B , respectively. Within network A, the nodes 
are randomly connected by A-cdgcs with degree distribu- 
tion Pa (A:), while the nodes in network B are randomly 
connected by B-edges with degree distribution Pg(fc). In 
our model a fraction qA of network A nodes depends on 
the nodes in network B and a fraction q B of network B 
nodes depends on the nodes in network A. We find that 
for strong coupling (large values of qA and q B ) the net- 
works undergo a first order transition while for a weak 
coupling they undergo a second order phase transition. 
Even for the case of weak coupling in which a second or- 
der percolation transition occurs, the system still disinte- 
grates in an iterative process of cascading failures unlike 
a regular second order percolation transition for a single 
network. 

The iterative process of cascading failures starts with 
randomly removing a fraction 1 — p of network A nodes 
and all the A-edges that are connected to them. Due to 
the interdependence between the networks, the nodes in 
network B that depend on removed A-nodes are also re- 
moved together with the B-edges that are connected to 
them. As nodes and edges are removed, each network 
breaks up into connected components, that we call clus- 
ters. We assume that when the network is fragmented, 
the nodes belonging to the giant component connecting 
a finite fraction of the network, are still functional, while 
nodes that are parts of the remaining small clusters be- 
come non-functional. Since each network is connected 
differently, the nodes that become non-functional on each 
step are different for both networks. This leads to the re- 
moval of more dependent nodes from the coupled network 
and so on. 

Next we present the formalism for the cascade process 
step by step. We define pa and pb as the fraction of 
nodes belonging to the giant components of network A 
and B respectively. The remaining fraction of network 
A nodes after an initial removal of 1 — p is a[ = p. The 
initial removal of nodes will disconnect additional nodes 
from the giant cluster. The remaining functional part of 
network A therefore contains a fraction a\ = a' 1 pA(o^' 1 ) 
of the network nodes. Since a fraction q B of nodes from 
network B depend on nodes from network A, the number 
of nodes in network B that become non functional is (1 — 
= <Zb(1 — a 'iP a{ol'i)) ■ Accordingly, the remaining 
fraction of network B \s j3[ = 1 — q B (l — a^pAipc'ij) and 
the fraction of nodes in the giant component of network 
Bis/?! =fip B (Pi). 

Following this approach we can construct the sequence, 
a n and /3 n , of giant components, and the sequence, a' n 
and f3' n , of the remaining fraction of nodes at each stage 
of the cascade of failures. The general form is given by: 
a[ =p, a\ = aipA(ai), 
P[ = 1 - q B (l - p A («i)p), ft = P'iPb(P'i), 
a' 2 = l- a[[l - q A (l - pb(P[))], «2 = a' 2 p A (a' 2 ) . . . 
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FIG. 2: An iterative process of failures for ER networks of 
size Na = Nb — 8 x 10 5 (a) A first order iterative process 
for p = 0.7455, a = b = 2.5, q A = 0.7 and q B = 0.6. (b) 
A second order iterative process for p — 0.605, a — b = 2.5, 
qa = 0.2 and qB = 0.75. Symbols represent simulation re- 
sults for different random realizations of the networks. Solid 
lines represent the solution of system @. (c) The fraction 
of nodes in network's B giant component, /3oo, as function 
of q A computed at p = pi(q A ), the line of the first order 
phase transition. The results are obtained by solving sys- 
tem (0) with additional condition dfA/dfs X dfs/dfA = 1 for 
a — 3, b — 3, qs = 1 • Inset: The same results (solid line) as 
function of \q A — q Ac \ yield a straight line with slope /3 = 1 
in double logarithmic scale. If qA is changed but p — p c is 
kept constant we obtain a straight line with slope /3 = 0.5 
(dashed line), (d) Simulation results for the phase transition 
of /Soo as a function of p for TV = 50K. For strong coupling 
between the networks we observe a jump in floo as expected 
in the first order phase transition (ER(circle) and SF(up rect- 
angle)). For weak coupling between the networks the change 
in /3oo is gradual as expected for the second and higher order 
phase transitions (ER(square) and SF(right rectangle)). 

a ' m =P\ l - Qa(1 ~ PB{fi' m ))], a m = a' m p A (a' m ), 

P'm = 1 - QB(l-PA(a' m )p), j3 m = P'nPBiP'm)- 

To determine the state of the system at the end of the 
cascade process we look at /3' m and a' m at the limit of m — > 
oo. This limit must satisfy the equation a' m —a' m+1 (or 
P'm—P'm+i) since eventually the clusters stop fragmenting 
and the fractions of randomly removed nodes at step m 
and m + 1 are equal. Denoting j3' m = y and a' m — x we 
arrive to a system of two equations with two unknowns: 

f y = i- 3s (1 ~Pa{x)p) m 
\x=p[l-q A (l-p B (y))}. {) 

The model can be solved analytically using the appara- 
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tus of generating functions. The generating functions 
will be denned for network A while similar equations 
describe network B. As in Refs. 11, Il3 we will intro- 



duce the generating function of the degree distributions 
Gao(0 = Sfc-^4(^)£ fe - Analogously we will introduce 
the generating function of the underlining branching pro- 
cesses, Gai(£) = G' AQ {£/) / G' AQ {1) . Random removal of 
fraction 1 — p of nodes will change the degree distribu- 
tion of the remaining nodes, so the generating function of 
the new distribution is equal to the generating function 
of the original distribution with the argument equal to 
1 — p(l — £) [ll|. The fraction of nodes that belong to 
the giant component after the removal of 1 — p nodes is 

p A (p) = l-G A0 [l-p(l-f A )}, (2) 

where f A — Ia{p) satisfies a transcendental equation 

f A = G A1 [l-p(l-f A )}. (3) 

In case of two ER networks, whose degrees are Poisson- 
distributed [II [11], the problem can be solved explicitly. 
Suppose that the average degree of the network A is a and 
the average degree of the network B is b. Then, G A \ (£) = 
G A0 = exp[a(£ - 1)] and G B1 (0 = G B0 = exp[6(£ - 1)]. 
Accordingly, p B {x) = 1 — /b and p A (x) = 1 — f A and 
therefore system (1) becomes 



X = p[l - q A f B ] 
y=l-q B (l-p[l-f A }), 



(4) 



where f A and f B satisfy the transcendental equations 



f A = aq>[ax(f A - 1)] 
f B = cxp[by(f B - 1)]. 



(5) 



The fraction of nodes in the giant components of net- 
works A and B respectively, at the end of the cascade 
process are given by a M = p(l - f A )(l - q A f B ) and 
Poo = (l-.f B )(l + q B (l-p)-pq B f A ). Fig.© shows ex- 
cellent agreement between computer simulations of the 
cascade failures and the numerical results obtained by 
solving systems ((4]) and ([5]). Excluding x and y from 
systems Q and ([5]), we obtain a system: 



f A = e -ap{fA-l)(qAfB-l) 

f B = e -b(q B (l-p[l-fA])-l)(fB-l)_ 



(6) 



The first equation can be solved with respect to f B and 
the second equation can be solved with respect to f A 



J B q A 



1 



h^-T^rjlfA^i; V/b, f A 



bp(f B 



: 1 

(7) 

The solutions of system ([7]) can be graphically pre- 
sented on a f A ,f B plane (Fig. [3]). The solutions are 
presented as a crossing of either f B (f A ) or f A = 1 
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FIG. 3: Illustrations of the different graphical solutions of 
system (O (see manuscript for detailed explanation of the 
different plots). 



with f B (f A ) or f A — 1 and are restricted to the square 
0</a<1;0</ b <1. There are three different 
possible solutions: (i) The solution where the giant com- 
ponents of both networks are zero (f A =l and f B =l) as 
in Fig|3](c). (ii) A solution for which only one of the gi- 
ant components of either network A or B is zero ( f A = 1 
and f B ^ l or f A ^ l and f B = l) as in FiglHd) (or 
FigJ3Je)). (in) A solution for which both networks have 
a non-zero giant component (f A ^ 1 and f B =^ l). This 
solution is given by the lowest intersection point of the 
curves in Figl3][a). This solution may disappear in two 
different scenarios. 

The first scenario is presented in Figl3^b) in which an 
infinitesimal change Az in the vector of the system pa- 
rameters z = {a,b,q A ,q B ,p) may lead to a first order 
phase transition in which the size of one or both of the gi- 
ant components changes discontinualy from a finite value 
to zero: (Figl^a) ->■ FigE^b) FigO|c) or Figl^d), or 
FigEfe)). The condition for the first order phase tran- 
sition is d ^f/^ rf ^/f" B ' > = 1 corresponds to the touching 
point of the two curves as in FigJ3Jb). When adding this 
condition to the two equations in system ([7]) we can find 
the three unknowns f A — f Al ,f B = f Bl and p = pi for 
given a,b,q A ,q B . Fixing a,b,q B will define a first or- 
der phase transition line p = pi(q A ) as function of q A 
[Fig. Kb)]. 

The second scenario is presented in Fig|3jf). In this 
case (corresponding to f A < 1, f B — 1 or equivalently to 
q B > 1 — 1/6), flao, continually decreases to zero, while 
a x stays finite. This situation corresponds to the second 
order phase transition that can be found by substituting 
f B = 1 into system ([7]). These two equations allow one to 
find f A = f AlI , and p = pjj which for fixed a, b, q B define 
a line of second order phase transitions p = pu(q A ) as a 
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FIG. 4: The critical point parameters, p c and qA„, as func- 
tions of the coupling strength qs for ER networks are plotted 
for different values of a = b. (a) For qs = 1 the networks 
(with large degrees) have the same critical coupling strength 
qA c = W[exp(3)] = 0.20794. The networks with small degrees 
do not have first order phase transitions (no critical points for 
qs), because for large values of qs, p c (gs) > 1, which is un- 
physical. The range of qs values for which the first order 
phase transition exists shrinks as a — b decreases and eventu- 
ally disappears for a = b — 3/2, when the critical point exists 
only for qA c = qs = 1/3 and p c = 1. This point is marked by 
a solid circle. 



changes only qA, then lim x ^ (l — /b)A/# = C" < 
corresponding to j3 — 1/2. The inset of Fig.[2jc) confirms 
our analytical predictions numerically. 

Although our analytical theory is developed for ER 
networks, the same qualitative conclusions hold for ran- 
domly connected networks with arbitrary degree distri- 
butions, since functions pa{x) and pb(u) can be ex- 
pressed in terms of generating functions of these distri- 
butions. Hence an analysis similar to Fig. 3 holds for any 
degree distributions. Computer simulations of interact- 
ing SF networks and ER networks presented in Fig. [2jd) 
support this analysis. 

We thanks the European EPIWORK project, the Israel 
Science Foundation and the DTRA for financial sup- 
port. S.V.B. thanks the Office of the Academic Affairs 
of Yeshiva University for funding the Yeshiva University 
high-performance computer cluster and acknowledges the 
partial support of this research through the Dr. Bernard 
W. Gamson Computational Science Center at Yeshiva 
College. 



function of qa [Fig.QJb)]. 

The line of the first order phase transitions merges with 
the line of the second order phase transitions in a crit- 
ical point which can be found by adding to system ([7]) 
both the first order condition dfB }/ A ^ d MM = \ an d 

djA d.JB 

the second order condition /g = 1 or /a = 1. These 
four equations allow us to find the critical parameters 
!b = !b c or f A = f Ac ,P=Pc and q A = qA c as functions 
of a,b,qB- Fig. [4] presents the solution for pdqs) and 
Qa c (qb) for different values of a(= b). The kink in the 
solutions occurs when both curves tangentially intersect 
at /a = 1)/b = 1 which corresponds to q~s = 1 — 1/b. 
The minimal value of p c occurs exactly at the kink, defin- 
ing the condition for the first order phase transition as 
Pc(q~b) < 1- Thus the first order transition can exist only 
in dense networks with sufficiently high average degrees, 
such that 4 (a — — 1) > 1. Low degree networks must 
disintegrate in the second order phase transitions. 

At the critical point the system can be reduced to a 
single transcendental Lambert equation. For the most 
simple case a/b = qB = 1, we find that Ja c — l/z> qA c — 
z — 2, p c = z/[a(z — 1)] and ceoo = (3 — z)/a, where 
z = W[exp(3)] — 2.20794 satisfies the Lambert equation 
z exp(z) = exp(3). 

To find the critical exponent /? near the critical point 
we express the order parameter Poo(<1a) as function of 
1A > q_A c along the transition line p = pi(qA) (inset of 
Fig. El[c)). Expanding fs in series of x = qA — qA c we 
find that lim x _ > o(l — /b)/x — C > 0, indicating that 
(3 = 1. Interestingly, if one keeps p = p c constant and 
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